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Abstract. The two-dimensional ferromagnetic Kondo model with classical 
corespins is studied via unbiased Monte Carlo simulations for a hole doping up 
to x = 12.5%. A canonical algorithm for finite temperatures is developed. We 
show that with realistic parameters for the manganites and at low temperatures, 
the double-exchange mechanism does not lead to phase separation on a two- 
dimensional lattice but rather stabilises individual ferromagnetic polarons for this 
doping range. A detailed analysis of unbiased Monte Carlo results reveals that 
the polarons can be treated as independent particles for these hole concentrations. 
It is found that a simple polaron model describes the physics of the ferromagnetic 
Kondo model amazingly well. The ferromagnetic polaron picture provides an 
obvious explanation for the pseudogap in the one- particle spectral function A^lw) 
observed in the ferromagnetic Kondo model. 
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1. Introduction 

Manganese oxides such as Lai-^Sr^MnOa, Lai-^Ca^MnOa and Laa—asSri-^ajMiiaOV, 
which have been thoroughly studied due to their colossal magnetoresistance (CMR), 
show a very rich phase diagram depending on doping, temperature, pressure and 
other parameters, see e.g [1,2]. The phase diagram includes ferromagnetic (FM), 
antiferromagnetic (AFM), paramagnetic (PM), charge-ordered, metallic as well as 
insulating domains. The manganites crystallise in the perovskite like lattice structure, 
and quasi two-dimensional systems with well separated Mn02-(bi)layers also exist. 

Crystal field splitting divides the five d-orbitals of the Mn ions into three 
energetically favoured ti g and two e g orbitals. All three ti g orbitals are singly occupied 
and rather localised. The filling of the e g orbitals is determined by doping and 
these electrons can hop from one Mn ion to the next via the intermediate oxygen 
(double exchange, DE). Due to a strong Hund's rule coupling, the spins of the 
three ti g electrons are aligned in parallel and form a corespin with length S = 3/2. 
Beinig localised, these electrons interact through superexchange which leads to a weak 
antiferromagnetic coupling between the corespins. Hund's rule coupling also leads to 
a ferromagnetic interaction between the itinerant e g electrons and the t2 g corespin. 
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Apart from double- and superexchange, a complete description would have to 
include Coulomb repulsion between the e g electrons, lattice distortions and disorder. 
Full quantum mechanical many-body calculations for a realistic model, including all 
degrees of freedom, are not yet possible, see however [3] for a one-dimensional study. 
Several approximate studies of simplified models have therefore been performed in 
order to unravel individual pieces of the rich phase diagram of the manganites. The 
electronic degrees of freedom are generally treated by a Kondo lattice model [4]. 

As full quantum mechanical results in more than one dimensions are difficult to 
obtain for the quantum mechanical Kondo lattice model, it has been proposed to treat 
the S=3/2 corespins classically, see de Gennes [5], Dagotto et al. [6,7] and Furukawa [8] 
Unbiased Monte Carlo (MC) techniques can then be applied. Coulomb repulsion [9], 
classical phonons [10] and disorder [11] have also been treated within this classical 
approach. The validity of this approximation has been tested in [6, 12-14] and it 
appears that quantum effects are important for (S=l/2) corespins or at T = 0. For 
finite temperature and S=3/2, classical spins present a reasonable approximation. 

Further approximations can be made by taking into account, that the Hund 
coupling Jh is much stronger than the kinetic energy. Consequently, configurations 
are very unlikely in which the electronic spin is antiparallel to the local corespin. A 
customary approach is to take Jh — > oo. This approximation however breaks down 
for the almost completely filled lower Kondo band. In the dilute hole regime, the full 
Kondo model is governed by an effective AFM interaction between the corespins due 
to excitations into the upper Kondo band. This effect is completely absent from the 
Jh — > oo model. 

An effective spinless fcrmion (ESF) model [15] has been proposed to improve 
upon this approximation. In this model, virtual excitations account for effects of 
configurations, in which the itinerant electron spin is antiparallel to the local corespin. 
It has been demonstrated that the results of the ESF model are in excellent agreement 
with those of the original Kondo model even for moderate values of Jh • 

For the FM Kondo with classical t2 g corespins, elaborate Monte Carlo simulations 
have been performed in various dimensions [6-8, 15-23] in order to determine the 
physical properties of the DE model. For a review see for example [7] and references 
therein. A two-dimensional Kondo lattice model for manganites in the Jh — > oo limit 
has been thoroughly investigated [23] by means of MC calculations similar to ours 
and by analytical comparison of the groundstate energy for several phases. Using a 
relatively high value for the antiferromagnetic exchange coupling, Aliaga et al. find 
phase separation (PS), stripes, island phases (small ferromagnetic domains that are 
stacked antiferromagnetically) for commensurate fillings, and a so called "Flux Phase". 
In a 2D-Kondo model applied to cuprates, stripes and a pseudogap are observed in 
MC simulations [24-26]. 

Many of these studies revealed features, e. g. an infinite compressibility near the 
filled lower Kondo band, which have been interpreted as signatures of PS. PS has 
also been reported [27] from computations based on a dynamical mean field treatment 
based on the DE model at T = 0. In previous MC studies [15, 21] for the DE model 
with classical core spins for ID systems, we had obtained numerical data comparable 
to those reported in Refs. [6,16,17]. A detailed analysis of the data [22] revealed 
however, that the aforementioned model with the standard parameter set, relevant for 
the manganites, favours individual polarons over phase separation. Other authors 
also found ferromagnetic polarons for the almost empty lower Kondo band (i. e. 
very few electrons) for S=l/2 corespins [28,29], for the ID AF Kondo Model with 
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few electrons [30], and for the ID paramagnet at higher temperatures [31]. Small 
ferromagnetic droplets were predicted from energy considerations [32]. 

In this paper, we present a numerical study of the 2D ferromagnetic Kondo model 
with classical corespins. As in ID, we find that the correct physical interpretation 
of the features which have been interpreted as PS is rather given by ferromagnetic 
polarons, i.e. small FM regions with one single trapped charge-carrier. The polaron 
picture allows also a straight forward and obvious explanation of the pseudogap, 
which has been previously observed in the spectral density in experiments [33-36] 
and MC simulations [7,15]. Experimental evidence for small FM droplets in low 
doped Lai-zCa^MnOs has been reported [37,38]. 

This paper is organised as follows. In Sec. El the model Hamiltonian is presented 
and particularities of the MC simulation for the present model are outlined. A 
canonical algorithm is introduced. In Sec. we introduce a simplified model of 
ferromagnetic polarons embedded in an AFM background and present results for this 
model. In Sec. these results are compared to unbiased Monte Carlo simulations 
of the 2D ferromagnetic Kondo model with classical corespins at realistic parameter 
values. The key results of the paper are summarised in Sec. 

2. Model Hamiltonian and unbiased Monte Carlo 

In this paper, we will concentrate solely on properties of the itinerant e g electrons 
interacting with the local ti g corespins. We also neglect the degeneracy of the e g 
orbitals. The degrees of freedom of the e g electrons are then described by a single- 
orbital Kondo lattice model [21]. As the corespins are approximated by classical spins, 
they are replaced by unit vectors S^, parameterised by polar and azimuthal angles 0i 
and (pi, respectively. The magnitude of both corespins and e g -spins is absorbed into 
the exchange couplings. 

2.1. Effective Spinless Fermions (ESF) 

By choosing the quantisation of the e g spin parallel to the local ti g corespin, simplified 
low-energy models for fillings < N c \ < 1 (i.e. for the lower Kondo band) can be 
derived, namely the Jh —> oo approximation and the effective spinless fermion model 
with finite Jh (see [15]): 

h = - E 'S 4 s- - E %f 4c, + j' E * • s, . (i) 

<i,j> id <id> 

The spinless fermion operators c A correspond to local spin-up electrons (i.e. parallel 
to the corespin) only. The spin index has therefore been omitted. With respect to a 
global spin-quantisation axis the ESF model Q still contains contributions from both 
spin-up and spin-down electrons. 

The first term in Eq. Q corresponds to the kinetic energy in tight-binding 
approximation. The modified hopping integrals tl'J depend upon the ti g corespin 
orientation 

*id = to u id ' ( 2 ) 
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where the relative orientation of the ti g corespins at site i and j, expressed by the 
angles < d < ir and < 4> < 2ir, enters via 

vff(S)= c iCj + s lSj e^-^ = cos(i%/2) 
u^;'r a (S) = Cisp-^i + SiCje-W* = sin(%/2) e ix « U 

with the abbreviations Cj = cos(^i/2) and Si = sin(i?j/2). These factors depend on 
the relative angle of corespins S.j and S j and on some complex phases ip^j and \ij . 
Because hopping is largest for parallel corespins, this term favours ferromagnetism. 
For certain spin structures, an electron may obtain a different phase depending on the 
path taken from one lattice site to another. An example of such structures is the so 
called Flux phase [23,39,40]. 

The second term in Eq. accounts for virtual hopping processes to antiparallel 
spin-corespin configurations and vanishes in the limit Jh — > oo. For finite Jh, the 
ESF model takes into account virtual hopping processes to antiparallel spin-corespin 
configurations much in the same way as the t J-model includes virtual hopping to 
doubly occupied sites for the Hubbard model. Because this term is proportional to 
the density, it is most relevant near the completely filled lower Kondo band, where 
the kinetic energy is moreover reduced. The last term is a small antiferromagnetic 
exchange of the corespins. 

The hopping strength to will serve as our unit of energy. Jh is usually taken to 
be of the order of magnitude of 4ig to 8io, and J' of the order of to/100. 



2.2. Grand Canonical Treatment 

We define the grand canonical partition function as 

Z(ji) = JV[S] trce-^ 5 )-^) 

where tr c indicates the trace over fermionic degrees of freedom at inverse temperature 
[3 and chemical potential /i, and N is the operator for the total number of e g electrons. 
As H is a one-particle Hamiltonian, the fermionic trace can easily be carried out using 
the free fermion formula, yielding the statistical weight of a corespin configuration S 

= — " (5) 

Equation 10} is the starting point of grand canonical Monte Carlo simulations of 
the Kondo model [6] where the sum over the classical spins is performed via Markov 
chain importance sampling; the spin configurations S are sampled with the probability 
determined by the weight factor w(S\(i). 

In order to avoid the CPU expensive full diagonalization of the one-particle 
Hamiltonian, Motome and Furukawa [41] suggested to replace it with an expansion 
in Chebychev polynomials. The CPU time then scales with the system size L as 
0(L 2 log(L)) instead of 0(L 3 ) as for the full diagonalization and the algorithm 
can be easily parallelised. We found, however, that on single processors the full 
diagonalization can be accelerated to be faster than this approach up to system sizes 
of 10 3 lattice sites, while both algorithms would be too slow for the study of larger 
systems on present day's processors. The key for a faster full diagonalization is to 
exploit the structure of the Hamiltonian: The lattice sites are relabeled in order to 
obtain a band matrix with as few diagonals as possible. This is only an alternative 
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assignment of the linear index to the two dimensional lattice vector and does therefore 
not introduce any approximation or error. Fast library routines for band matrices can 
then be used. Since the purpose of our MC code was to perform parameter studies, 
we did not program a parallel algorithm but instead ran the whole calculation on 
each CPU with a different parameter set. Recently, an O(L) algorithm has been 
proposed by the same authors, which reduces the numerical effort by approximating 
the matrix- vector multiplication [20]. 

In the 2D case we have employed MC updates in which single spins were rotated. 
The angle of rotation was optimised to keep the acceptance high enough. From time 
to time a complete flip Sj — > — Sj was proposed. The skip between subsequent 
measurements was chosen to be 50 to a few hundreds of lattice sweeps reducing 
autocorrelations to a negligible level. We have performed MC runs with some 
hundreds to 2000 measurements on a 12 x 14 lattice. This geometry was chosen 
to reduce finite size (closed shell) effects observed on a square lattice. The number of 
measurements was higher for calculations in the polaronic regime, where the particle 
number fluctuates strongly in order to have sufficient measurements for each filling. 

As previously shown [21], the spin-integrated one-particle Green's function in 
global quantisation can be written as 



E 



« a ia - at CT » w = / V[S] w(S\v)u]} (S) « c. ( ; c\ »£ , (6) 



where <C c i\ c \ IS the Green's function in local spin quantisation. It can 

be expressed in terms of the one-particle eigenvalues and the corresponding 
eigenvectors of the Hamiltonian H{S): 

^ 1 ' j >>UJ ^ u- (eW -n) + i0+ 

It should be pointed out that the one-particle density of states (DOS) is identical in 
global and local quantisation; for details see [21]. 



2.3. Canonical Algorithm 

Due to the jump in the electron density at the critical chemical potential (shown 
later in Fig. EJ) , some electron fillings cannot be examined with the grand canonical 
algorithm. We therefore developed a canonical scheme. Canonical calculations were 
done by computing the eigenenergies for each corespin configuration and then filling 
the available electrons into the lowest levels [23]. This method does, however, not 
account for thermal particle-hole excitations around the Fermi energy. When several 
one-particle states have similar energy, these may become important. It was also 
proposed to adjust the chemical potential by solving an implicit equation using the 
Newton-Raphson algorithm for every single spin configuration to give the wanted 
particle number [42]. While this approach includes particle hole excitations, it is still 
not certain, whether the calculation is correct and small differences possibly have a 
considerable effect when phases with and without a pseudogap compete. 

On the other hand, an exact approach would mean calculating the Boltzmann 
weight for every possible distribution of N e i particles on L energy levels and summing 
over their contributions. Even for small lattice sizes L, this clearly becomes too 
demanding for more than a few electrons or holes. Instead, we took into account just 
the lowest excitations of the Fermi sea by filling < N e i electrons into the lowest 
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states and considering only the distributions of the N e i — N®, remaining electrons on 
the states around the Fermi energy. Usually, it is sufficient to take N e [ — N®, ~ 5. The 
weight for the corespin configuration S then depends on the particle number instead 
of the chemical potential: 

y . e -PH{S,V(N el )) 

'"^ = i ^-mr ) — • < 7 > 

where V denotes these restricted permutations. 

Although this is more time-consuming than the grand canonical calculation of the 
fermionic weight, the additional consumption of computer time is small compared to 
the time needed for the diagonalization of the one-particle Hamiltonian. The particle- 
hole excitations, which are thus included, can be crucial when examining competition 
between phases with and without a pseudogap. 

A MC update - especially a complete spin flip - may lead to a configuration 
which is very unlikely to occur at the given particle number, although it may be a 
good configuration for a different filling. A later MC move might then lead back to 
the original particle number, and these moves should improve autocorrelation. We 
therefore allow density fluctuations within a set of four to five particle numbers. In 
order to spend a comparable number of MC steps at each filling, prior weight factors 
g(N e i) were introduced and adjusted in a prerun, giving 

w{S)= £; w(S\N el ) g(N e i) . (8) 

The sum is taken over the set of allowed particle numbers N m i n < N e i < N max and 
w(S\N e i) is calculated according to Eq. Q. 

When evaluating observables for fixed electron number, one has to calculate the 
expectation value 

{ u )n b1 = — ^ — fQlM s — , (yj 
22s w ( s \ N ei) 



which can be rewritten as 



EsO(S) Nel ^§§f w(S) 



(o) Nel = — " -r' — — , (io) 

l^S w(S) W ^> 

Configurations S occur in the Markov chain with probability proportional to w(S); 
when the sum is taken over the configurations produced by the MC run, the 
expectation value therefore becomes 

{V) Ncl ,MC= — w(S\N cl ) ■ 

1^S,MC w(S) 



3. Ferromagnetic Polaron Model 

Near half filling of a single e g band, a tendency toward phase separation has been 
reported in various computational studies [6, 16, 17, 27]. 

In most cases, the existence of phase separation is inferred from a discontinuity 
of the electron density as a function of the chemical potential. At the critical chemical 
potential where this discontinuity is found, it is claimed that the system separates 
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(a) 



(b) 



Figure 1. Idealised FM polaron of Lf = 5 lattice sites, embedded in an AFM 
background, (a) Spin and hole-density configuration for the ground state. Empty 
(filled) squares represent spin down (up). Height represents hole density, (b) 
Contribution of the polaron to the one-particle spectral function. For visibility, 
the (5-peaks in the spectral density have been broadened to a width of 0.2. 



into FM domains of high carrier concentration and AFM domains of low carrier 
concentration. 

We have already shown that for a ID system [22] this picture is in general 
incorrect. In fact, what happens is that each single hole is dressed by a ferromagnetic 
cloud in which it delocalizes. The system can be well described by free quasiparticles 
consisting of a single hole plus a local (three to four-sites for ID, 5 sites for 2D) 
ferromagnetic well embedded in an AF background. Each of these added quasiparticles 
gains the same energy, which is exactly balanced by the energy to be paid for the 
critical chemical potential ji* . Hence the discontinuity of the particle number at low 
temperatures. 

Here we show that ferromagnetic polarons, i.e. single charge carriers surrounded 
by small ferromagnetic spin-clouds, are indeed formed when holes are doped into a 
completely filled lower Kondo band in 2D. In this section we discuss the properties 
of idealised two dimensional model polarons, whereas in the next section we compare 
our polaron-model to unbiased Monte Carlo results. 



3.1. One single polaron 

As reference configuration we consider the completely filled lower Kondo band where 
super exchange and the contribution from the virtual hopping process, see Eq. (Q, 
give rise to an antiferromagnetic corespin pattern. The smallest defect in a completely 
antifcrromagnetically ordered 2D lattice is the flipping of one single spin. Then a 
five-site ferromagnetic region forms. When we introduce a hole into the system, it 
can delocalize in this region because of the double-exchange mechanism. As there 
is no double-exchange hopping between sites with perfectly antiferromagnetic spins, 
the hole is trapped inside the ferromagnetic region. In this simple model, the hole 
can hop from the central site of this region (site 1) where the spin has been flipped, 
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to all nearest neighbours (sites 2-5). 
situation reads 



The tight-binding Hamiltonian describing this 



hole 



/ o 
1 
1 
1 

V 1 



i \ 





o ) 



with the hopping between ferromagnetic sites ff in this model. Since the 
Hamiltonian is symmetric with respect to rotations by tt/2, the ground state 
should also show this property and can thus be found in the space spanned by 
(1, 0, 0, 0, 0), (0, 1/2,1/2, 1/2, 1/2), where the Hamiltonian reads 

2 
2 



-ffholc — —tf 



1(a) 



2tf. This energy 
The ground state 
Excited states 



The derealization energy of the hole is thus given by ehoic = 
can be gained as a hole delocalizes in the ferromagnetic domain 
1/V2 (1, 1/2, 1/2, 1/2, 1/2) has the hole density depicted in Fig ~ 
are found at e = +2 if and e = 0. The highest state at e = +2 if is given by 
l/y/2 (—1, 1/2, 1/2, 1/2, 1/2) and does also have s-symmetry, while the states at e = 
have p- and d-symmetry. They appear in the one-particle spectral function of the 
configuration shown in Fig 



1(b) 



To create such a ferromagnetic domain, however, four antiferromagnetic bonds 
have to be broken. This costs the energy of 2 x 4 J c g. Near the completely filled lower 
Kondo band x = 0, the total antiferromagnetic exchange coupling is approximately 
given by J e g — 1/(2 Jh) + J', sec [22]. The energy gained by adding one hole to the 
system thus reads 



e P oi — —2tf + 8 J c ff 



(12) 



When the chemical potential approaches — e po i from above, holes start to enter the 
system forming individual polarons. Therefore, the critical chemical potential is given 
by pi* = -e po i. 

Up to a certain concentration, these holes can be treated as free fermions which 
all have the same energy e po i. The energy may, however, depend on the temperature if 
tf or J e ff do. The more obvious temperature effect is the smearing of the discontinuity 
in the electron filling at higher temperatures, which results from the application of 
the Fermi-Dirac statistics to these quasiparticles. 

It is interesting to compare the energy of two such independent holes to a larger 
FM formation containing two holes. The two possibilities with the least broken AFM 
bonds are depicted in Fig. For both configurations, eight AFM bonds have to be 
broken, the same as for two polarons. For the formation in Fig. 2(a) the sum of the 
two lowest eigenenergies is however t\ +e 2 = —3.97 if, i. e. the kinetic energy gained is 
smaller than for two independent polarons, where it is 2ei 10 i c = —4 tf. For the second 
formation Fig. 2(b) the energy gain is even lower (ei + e 2 = —3.86 if). The larger FM 
structures with two holes therefore have higher energy than two separate holes and 
are energetically disfavoured, although the difference is not very large, especially for 
the structure depicted in Fig. 2(a) 
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Figure 3. Representative hole-density configuration for 10 (left, x fs 0.06) and 
20 (right, x m 0.12) holes in the simplified polaron model on a 12 X 14 lattice. The 
holes are located where one spin is flipped from the AFM reference configuration. 
The results correspond to j3 = 50, J' = 0.02, Jjj = 6. Grey shades are for better 
visibility. 




r X Mr 



Figure 4. Spectral density of the polaron model for 20 holes on a 12 X 14 lattice 
(x ~ 0.12), using the same parameter values as Fig. 131 
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Figure 5. Electron density y as a function of the chemical potential fj, of a 14 X 12 
lattice for /3 = 50, J' = 0.02 and Jjj = 6. Circles indicate MC results, the solid 
(dashed) line shows the free fermion results for the fitted (estimated) energy. 

3.2. Results for a small number of polarons 

In this subsection we push our polaron ideas further to treat the case of a small number 
of holes in an AFM background. To this end, we perform a simple simulation. We 
start with a perfect AFM reference configuration. We then add iVhoies holes by flipping 
JVhoies randomly chosen spins, excluding flipping a spin back. Finally, we add random 
deviations to each corespin in order to account for thermodynamic fluctuations. These 
fluctuations lead to a finite DE hopping amplitude in the AFM band and their size 
is therefore fitted to match the bandwidth observed in the MC results, see Sec. 0] 
We then diagonalise the resulting ESF-Hamiltonian and compute observables in the 
canonical ensemble with = A^ites - A ho ics = A si tcs - ^flipped spins as explained in 
Sec. 12.81 The observables are averaged over many such configurations. Typical hole- 
density configurations are depicted in Fig.Olfor 10 and 20 holes, i. e. 10 resp. 20 spins 
were flipped from the initial perfect AF configuration. In the formulae for the ESF- 
Hamiltonian and the observables, the parameters were set to (3 — 50, J' = 0.02 and 
Jh = 6. It should be emphasised that the probability distribution used in choosing the 
spins to be flipped is completely flat. It is only ensured that exactly Ahoies spins are 
flipped from the AFM reference configuration. The chosen parameter values therefore 
do not influence the obtained spin configurations. 

As an example for observables calculated in this pure polaron model, we show 
the one-particle spectral function for 20 holes in Fig. 0] The center is occupied by 
the AFM tight-binding band with a mirror band due to the doubling of the unit cell 
in the AFM spin configuration. At u> — ±2, the polaronic states can be clearly seen. 
In addition to them, one sees a number of weaker signals in the vicinity of uj = ±2. 
These stem from larger ferromagnetic regions, i. e. from contiguous polarons. 

4. Unbiased MC Results 

In this section we present results of unbiased Monte Carlo simulations for J' = 0.02 
and Jh = 6 and show that they correspond to independent ferromagnetic polarons. 

Figure shows the electron density as a function of the chemical potential at 
= 50. Depending on the value of the hopping parameter t$, this is in a range of 
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Figure 6. MC snapshot of the hole density for 10 (left) resp. 20 (right) holes in 
a 14 X 12 lattice at /3 = 50, J' = 0.02, J H = 6. 




energy energy 



(a) (b) 

Figure 7. One-particle DOS for few holes in a 14 X 12 lattice, J' = 0.02, J H = 6. 
(a): Temperature-dependence of the DOS in the case of one hole. The inset 
shows an enlargement of the polaronic peak, (b) Doping dependence of the DOS 
for p = 50 

50X — lOOif , i. e. relevant for experiments. There is a discontinuity in the density 
(infinite compressibility), but one observes that the electron number does not drop 
at once from 1 (AFM) to 0.8 (FM). Instead, it first decreases only slowly from the 
completely filled band, the slope of the curve then becomes gradually steeper until it is 
vertical. For a qualitative description of the MC results by the polaron model, we use 
tf = t and J c ff = 1/(2Jh) + J' which yields a polaron energy of e po i ~ —1.17. Using 
this value for the critical chemical potential we obtain the dashed line. Although it is 
shifted by some constant energy, it already correctly reflects the trend in the Monte 
Carlo data. Much better agreement can be found by fitting the polaron energy to 
the Monte Carlo data. In our case we obtain e po i ~ —1.14. The corresponding Fermi 
function is shown as the solid line. 

Figure H3 shows MC snapshots with 10 and 20 holes. The polarons can be clearly 
seen, 10 polarons for 10 holes and 19 polarons for 20 holes. Only the 20 th hole at the 
larger doping is delocalised. There is an obvious similarity to the idealised polaron 
model, see Fig. 

As in the one-dimensional case, the polarons induce separate states in the one- 
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particle density of states (DOS) depicted in Fig. [7] as described in Sec. [3] Figure 7(a) 



shows the DOS in the case of one hole in a 14 x 12 lattice at various temperatures. 
One observes a broad peak in the center and two polaronic peaks at ±2 that are 
separated by a pseudogap. The pseudogap observed in the Kondo model is thus a 
direct result of the FM polarons. The broad central peak is due to holes moving in 
the not quite perfectly antiferromagnetic background. As a result of the superexchange 
term in the Hamiltonian, it is centered around e = — z/(2Jn) with z = 4 in the 2D 
square lattice. When no hole is in the system, this peak also shows up and is then 
the only feature of the DOS. The width of the peak is mainly dominated by corespin 
fluctuations around the completely ordered state. Consequently the width decreases 
with decreasing temperature. This leads to a wider pseudogap at lower temperatures, 



as depicted in Fig. 7(a) 



The polaronic peak, on the other hand, remains largely unaffected by 
temperature. As can be seen in the inset of Fig. 7(a) its shape is virtually constant. 
Upon introducing more holes (see Fig. 7(b) I, the weight of the polaronic peaks 
increases whereas their position and the shape of the central band remain unaffected. 
The weight of the polaronic peaks corresponds to the number of holes. The shape of 
both the antiferromagnetic band and the polaronic peaks only begins to change when 
a large number of holes are added, so that the probability for overlapping polarons 
becomes considerable. With 5 sites per polaron, 20 polarons would take 100 sites, or 
60% of a 14 x 12 lattice, therefore connected polarons and disturbances of the AF 
background are to be expected. It is therfore remarkable that even with 20 holes, the 
polarons seem to remain largely independent. 




(a) 



(b) 



Figure 8. Spectral density for J' = 0.02, f3 = 50, J H = 6 on a 12 X 14 lattice, 
(a): six holes (x S3 0.035): Polaronic states in addition to the AFM band, (b): 20 
holes (x S3 12%) 



Figure[H]shows the spectral density for the ESF model with Jh = 6, J' = 0.02 and 
j3 = 50 for 6 and 20 holes. In addition to the central band, the polaronic states can 
be seen at energies slightly below ±2, in perfect agreement with the simple polaron 
model, see Fig. 1(b) The states at oj ~ are lost within the AFM band. The results 
for 20 holes (x ss 0.12, Fig. 8(b) I are slightly smeared compared to the results for the 
simple polaron model in Fig. |3J but the similarities are striking. 

In order to verify that the addition of holes leads primarily to more small polarons 
rather than to a growth of the existing ones, we use a dressed corespin correlation 
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r 2 

Figure 9. Dressed corespin correlation Eq. 1131 for J' = 0.02 from unbiased 
MC Data for 1 (o), 6 (x) and 20 (A) holes. Continuous lines are data for the 
simple Polaron model (see Sec. [3J: 1 Polaron (dotted), 6 (solid) and 20 Polarons 
(dashed). The inset shows the corespin correlation S(r) = ■ S-*,^. The 

MC Simulations were done for Jjj = 6, = 50 on a 12 X 14 lattice. 

function 

i 

Where ni is the hole density at site i related to the electron density via ni ~ 1 — n^, 

the sum over i is taken over all lattice sites. This dressed correlation measures the 
ferromagnetic regions around holes. Figure El shows the results for 1, 6 and 20 holes 
on 12 x 14 sites, which correspond to doping levels of x — 0.006, x — 0.036, and 
x = 0.12. The results are almost independent of the doping level x and, above 
all, the ferromagnetic region does not grow with increasing hole density. The data 
are compared to those obtained for the independent polaron model introduced in 
Sec. |3 and show very good agreement. The inset of Fig. [5] shows the usual corespin 
correlation which reveals the AFM background. The antiferromagnetism decreases 
with increasing hole concentration, both for the MC simulations and the idealised 
polaron model. However, it shrinks somewhat faster for the full ESF model. 

The results begin to deviate from the independent polaron results only for 
doping levels x > 13%, with a more homogenous phase setting in at x ~ 21%. In 
between, the polarons attract each other with a tendency to form larger FM clusters 
in the antiferromagnetic background, which eventually leads to phase separation 
(PS). The transition from polarons to PS is not well defined. On increasing hole 
concentration the polarons first coexist with larger hole-rich clusters which then grow 
and finally dominate. Details on this range of doping and on the influence of the AFM 
superexchange parameter J' are subject of a subsequent publication. 

5. Conclusions 

In this paper, the ferromagnetic Kondo (double-exchange) model in 2D has been 
analysed by unbiased finite temperature Monte-Carlo simulations. It has been found 
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that upon hole doping, small ferromagnetic regions appear around each individual 
hole while the rest of the lattice stays antiferromagnetically ordered. Each of the 
ferromagnetic regions contains one single hole. Therefore, the physics close to half 
filling is not governed by phase separation into larger FM and AF regions, as previously 
reported, but by single-hole ferromagnetic polarons moving in an antiferromagnetic 
background. 

The critical chemical potential //* at which holes start to enter the lower Kondo 
band can be found from simplified energy considerations (Eq. 1)12(1 ^. For fj, significantly 
above /i*, the band is completely filled and the corespins are antiferromagnetic. 
Around fi* , holes enter the e g -band, forming isolated FM domains in the shape 
depicted in Fig. |l(a)| each containing one single hole. This is corroborated by MC 
snapshots, the functional dependence of the electron density on the chemical potential, 
the spectral density and the dressed corespin correlation Eq. (fTTTf) . 

The discontinuity in the electron density vs. the chemical potential (i. e. infinite 
compressibility) is usually taken as evidence for PS. In the case of the Kondo model, 
this discontinuity is a consequence of a large (macroscopic) number of degenerate 
polaron states. When the chemical potential is close to the energy of these states, 
the number of holes (polarons) in the lattice strongly fluctuates. The weight of the 
polaron peak in the spectrum is directly linked to the number of holes (Fig. 7(b) I. In 
order to obtain numerical results at a fixed hole number, it was necessary to develop 
a canonical algorithm for our Monte Carlo simulations. 

Another consequence of the formation of single-hole FM polarons is the opening 
of a pseudogap. The small FM regions of the polarons contain only a few electronic 
states that are energetically well separated from each other. Moreover, the width of the 
antiferromagnetic band is much smaller than the difference between the highest and the 
lowest polaron states. Therefore, no states can be found for energies between the upper 
edge of the antiferromagnetic band and the highest state within the polaron. This gives 
rise to a pseudogap in the one-particle spectral function. The same arguments explain 
the appearance of a mirror gap well below the chemical potential. 

A pseudogap is indeed observed in experiments [33-35] and in MC simulations for 
the Kondo model [7,15]. Experiments at low doped Lai-^Ca^MnOs showed evidence 
of small FM droplets in an AFM background [37,38]. 

Our analysis yields compelling evidence against the PS scenario and in favour 
of FM polarons for small doping in 2D for realistic parameter values for manganites. 
A similar behaviour has been previously found for ID. Furthermore, the coupling to 
lattice degrees of freedom will additionally localise holes and inhibit the formation 
of a ferromagnetic phase (Jahn- Teller polarons) [3, 12,43]. It should be noted that, 
depending on the value of the hopping parameter tg, the temperature investigated 
in this paper (/? = 50) is in a range of 50K — lOOi^T, which is in agreement with 
temperatures in experiments. 
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